import vtk
%matplotlib inline
from dipy.reconst.dti import fractional_anisotropy, color_fa
from argparse import ArgumentParser
from scipy import ndimage
import os
import re
import numpy as np
import nibabel as nb
import sys
import matplotlib
# matplotlib.use('Agg') # very important above pyplot import
import matplotlib.pyplot as plt
from dipy.reconst.dti import from_lower_triangular
img = nb.load('fibers/dogsigma_0gausigma_0tensorfsl.nii')
data = img.get_data()
output = from_lower_triangular(data)
output_ds = output[:, :, :, :, :]
print output.shape
print output_ds.shape
print data.shape
# FA = fractional_anisotropy(output_ds)
# FA = np.clip(FA, 0, 1)
# FA[np.isnan(FA)] = 0
# print FA.shape
from dipy.reconst.dti import decompose_tensor
evalues, evectors = decompose_tensor(output_ds)
print evectors[..., 0, 0].shape
print evectors.shape[-2:]
print evalues.shape, evectors.shape
# print FA[:, :, :, 0].shape
print evectors.shape
FA = fractional_anisotropy(evalues)
print FA.shape
RGB = color_fa(FA, evectors)
nb.save(nb.Nifti1Image(np.array(255 * RGB, 'uint8'), img.get_affine()), 'tensor_rgb_upper.nii.gz')
def plot_rgb(im):
plt.rcParams.update({'axes.labelsize': 'x-large',
'axes.titlesize': 'x-large'})
if im.shape == (182, 218, 182):
x = [78, 90, 100]
y = [82, 107, 142]
z = [88, 103, 107]
else:
shap = im.shape
x = [int(shap[0]*0.35), int(shap[0]*0.51), int(shap[0]*0.65)]
y = [int(shap[1]*0.35), int(shap[1]*0.51), int(shap[1]*0.65)]
z = [int(shap[2]*0.35), int(shap[2]*0.51), int(shap[2]*0.65)]
coords = (x, y, z)
labs = ['Sagittal Slice (YZ fixed)',
'Coronal Slice (XZ fixed)',
'Axial Slice (XY fixed)']
var = ['X', 'Y', 'Z']
idx = 0
for i, coord in enumerate(coords):
for pos in coord:
idx += 1
ax = plt.subplot(3, 3, idx)
ax.set_title(var[i] + " = " + str(pos))
if i == 0:
image = ndimage.rotate(im[pos, :, :], 90)
elif i == 1:
image = ndimage.rotate(im[:, pos, :], 90)
else:
image = im[:, :, pos]
if idx % 3 == 1:
ax.set_ylabel(labs[i])
ax.yaxis.set_ticks([0, image.shape[0]/2, image.shape[0] - 1])
ax.xaxis.set_ticks([0, image.shape[1]/2, image.shape[1] - 1])
plt.imshow(image)
fig = plt.gcf()
fig.set_size_inches(12.5, 10.5, forward=True)
return fig
affine = img.get_affine()
fa = nb.Nifti1Image(np.array(255 * RGB, 'uint8'), affine)
im = fa.get_data()
# print np.asarray(fa)
fig = plot_rgb(im)
import os
from PIL import Image
im = plt.imread('fibers/v100/ch0/luke40.tiff')
plt.imshow(im)
import dipy.reconst.dti as dti
from dipy.reconst.dti import fractional_anisotropy
from dipy.data import default_sphere
from dipy.direction import DeterministicMaximumDirectionGetter
from dipy.io.trackvis import save_trk
from dipy.reconst.csdeconv import ConstrainedSphericalDeconvModel
from dipy.data import get_sphere
sphere = get_sphere('symmetric724')
from dipy.reconst.dti import quantize_evecs
peak_indices = quantize_evecs(evectors, sphere.vertices)
print FA.shape
from dipy.tracking.eudx import EuDX
eu = EuDX(FA.astype('f8'), peak_indices, seeds=50000, odf_vertices = sphere.vertices, a_low=0.2)
tensor_streamlines = [streamline for streamline in eu]
len(tensor_streamlines)
tensor_streamlines_trk = ((sl, None, None) for sl in tensor_streamlines)
from argparse import ArgumentParser
from dipy.viz import window, actor
import numpy as np
def visualize(fibers, outf=None):
"""
Takes fiber streamlines and visualizes them using DiPy
Required Arguments:
- fibers:
fiber streamlines in a list as returned by DiPy
Optional Arguments:
- save:
flag indicating whether or not you want the image saved
to disk after being displayed
"""
# Initialize renderer
renderer = window.Renderer()
# Add streamlines as a DiPy viz object
stream_actor = actor.line(fibers)
# Set camera orientation properties
# TODO: allow this as an argument
renderer.set_camera() # args are: position=(), focal_point=(), view_up=()
# Add streamlines to viz session
renderer.add(stream_actor)
# Display fibers
# TODO: allow size of window as an argument
window.show(renderer, size=(600, 600), reset_camera=False)
# Saves file, if you're into that sort of thing...
if outf is not None:
window.record(renderer, out_path=outf, size=(600, 600))
visualize(tensor_streamlines)
from IPython.display import Image
def vtk_show(renderer, width=400, height=300):
"""
Takes vtkRenderer instance and returns an IPython Image with the rendering.
"""
renderWindow = vtk.vtkRenderWindow()
renderWindow.SetOffScreenRendering(1)
renderWindow.AddRenderer(renderer)
renderWindow.SetSize(width, height)
renderWindow.Render()
windowToImageFilter = vtk.vtkWindowToImageFilter()
windowToImageFilter.SetInput(renderWindow)
windowToImageFilter.Update()
writer = vtk.vtkPNGWriter()
writer.SetWriteToMemory(1)
writer.SetInputConnection(windowToImageFilter.GetOutputPort())
writer.Write()
data = str(buffer(writer.GetResult()))
return Image(data)
renderer = window.Renderer()
# Add streamlines as a DiPy viz object
stream_actor = actor.line(tensor_streamlines)
# Set camera orientation properties
# TODO: allow this as an argument
renderer.set_camera() # args are: position=(), focal_point=(), view_up=()
# Add streamlines to viz session
renderer.add(stream_actor)
vtk_show(renderer)
import IPython.display import Image
Image("Screenshots/vtk_face1.PNG", width=200, height=200)
Image("Screenshots/vtk_face2.PNG", width=200, height=200)
Image("Screenshots/vtk_face3.PNG", width=200, height=200)
Image("Screenshots/vtk_face4.PNG", width=200, height=200)
Image("Screenshots/vtk_top.PNG", width=200, height=200)
Image("Screenshots/vtk_bottom.PNG", width=200, height=200)
len(tensor_streamlines)
# tensor_streamlines[0][:,0:2]
test = tensor_streamlines[0:1000]
print len(test)
plt.figure(1)
plt.subplots(figsize=(10, 10))
plt.subplot(311)
plt.title("Y-axis vs X-axis (1000 fibers)")
for i in range(len(test)):
plt.plot(test[i][:,0], test[i][:,1])
plt.subplot(312)
plt.title("Z-axis vs X-axis (1000 fibers)")
for i in range(len(test)):
plt.plot(test[i][:,0], test[i][:,2])
plt.subplot(313)
plt.title("Z-axis vs Y-axis (1000 fibers)")
for i in range(len(test)):
plt.plot(test[i][:,1], test[i][:,2])
plt.tight_layout()
plt.show()
print len(tensor_streamlines)